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I. INTRODUCTION 

Multifractal models were first introduced in the 1960s 
by the so-called "Russian school" in turbulence theory 
[TJI2]. In turbulence, multifractality can be conceived as a 
weakening of the spatial selfsimilarity in the velocity field 
implicitly assumed in Kolmogorov's 1941-theory [3]. This 
generalization is called the Kolmogorov-Obukhov model 
and entails modeling the spatial variability of the energy 
dissipation rate as a random measure with certain multi- 
scaling properties. The Kolmogorov-Obukhov model is 
treated rigorously by Kahane and this construction is 
known as Gaussian multiplicative chaos. 

In recent years multifractal random processes and mul- 
tifractal random measures have received increased at- 
tention and are widely used in physics, geophysics and 
complex systems theory. Examples include phenomena 
as diverse as internet traffic [S], geomagnetic activity 
[6l [7] and rainfall patterns [8]. In addition, multifrac- 
tal processes provide natural models for the long-range 
volatility persistence observed in financial time series. 
This was first discovered by Ghashghaie [H] and Man- 
delbrot [1^, and since the late 1990s much work has 
been done on multifractal modeling of financial markets 
[TTl [T^ . Logarithmic returns of assets are modeled as 
Xt = X{t + At) — X{t), where X{t) are continuous- 
time processes with stationary increments and multifrac- 
tal scaling. The latter means that the moments of X{t) 
are power-laws as functions of time; 

E|X(<)|« (1) 

either in some interval t e (0, R) or asymptotically as 
i — >■ 0. The scaling function ({q) is linear for self-similar 
processes, but may in general be concave. Processes sat- 
isfying equation ([I]) with strictly concave scaling func- 
tions are generally referred to as multifractal. 

Two well-known "stylized facts" of financial time series 
are that log-returns are uncorrelated and non-Gaussian. 
Based on this, Mandelbrot \i3l deduced that if prices are 
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described as selfsimilar processes, then these processes 
must be so-called Levy flights, i.e. a-stable Levy pro- 
cesses with a < 2. However, if one allows non-linear 
scaling functions, then one can maintain uncorrelated 
log- returns by simply imposing the condition C(2) — 1. 
The concave shape of (^{q) implies that the variables X{t) 
are increasingly leptokurtic with decreasing t, and conse- 
quently non-Gaussian. Moreover, as opposed to Levy 
flights, multifractal processes have strongly dependent 
increments and can therefore describe a third "stylized 
fact" of financial time series, namely volatility cluster- 
ing. 

Notwithstanding that multifractal processes provide 
accurate and parsimonious descriptions of temporal fi- 
nancial fluctuations, the models are rarely implemented 
for forecasting and risk-analysis in financial institutions. 
This is partially due to a lack of accurate, stable and ef- 
ficient inference methods for multifractal processes. Pa- 
rameter estimation has so far mostly been made using 
various moment-based estimators, such as the general- 
ized method of moments (GMM). Alternatively, one can 
fit the estimated scaling functions to theoretical expres- 
sions of Cil)- However, as pointed out in e.g. jl4| and 
|15) . the standard estimators of scaling exponents have 
large mean square errors for time series of length compa- 
rable to those typically available in econometrics. 

An exception to the statements above is the Markov- 
Switching Multifractal (MSM) model [H] where maxi- 
mum likelihood estimation is feasible. In discrete time 
MSM implies that the increments Xt are described by a 
stochastic volatility model on the form 

Xt^c^^/Mtet. (2) 

Here et ^ A/'(0, 1) are independent variables and the 
volatility is a product on the form 

Mt ^ Mt,iMt,2 ■ ■ ■ Mt,K , 

where (for each time step t) Mt^k are updated from a 
distribution M with a probability 7^. = 1 — (1 — 71)^'" \ 
In this model however, maximum likelihood estimation is 
only possible in the case where M is defined on a discrete 
state space, and there is a limitation on the magnitude of 
K which should not exceed « 10 [T^. These restrictions 
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not only limit flexibility with respect to the distribution 
of returns, but also the possible range in the volatility 
dependency. 

This paper concerns parametric inference for the mul- 
tifractal random walk (MRW) introduced by Bacry et 
al. [17]. The increment process xt ^ X(t + At) — X{t) is 
still a discrete-time stochastic process described by equa- 
tion ([2]) , but now the volatility is modeled as Mt = c e'*' , 
where ht is a stationary and centered Gaussian process 
with the co-variance structure 

Cov( ht, h,) ^ X^\og+ f- — — , (3) 

^ ^ ^ {\t-s\ + l)At ^ ' 

where log^ a max{loga,0}. Here T > is called the 
correlation range [15] and A > is called the intermit- 
tency parameter. The constant c ensures normalization 
and is chosen so that 1/c = E[e''*]. We denote R = T/At. 

Let = (A, cr, R) denote the parameter vector and y = 
(2/1, ... , Un) € K" a fixed time series. The main result of 
this paper is the development of a method for efficiently 
computing approximations to the likelihood function 

Ci9\y)^pAy\9), 

where Px{'\9) is the probability density function of a 
random vector x = {xi, . . . ,Xn) produced by the MRW 
model with parameters 9. Using the likelihood function, 
parameters can be determined by means of the maximum 
likelihood (ML) estimator: 

9 — argmaxg £{9\y) . 

Our method exploits that the discrete MRW model has 
a construction similar to simple volatility (SV) models. 
The distinguishing feature is that the processes ht are au- 
toregressive in SV models. By truncating the dependency 
structure in the logarithmic volatility ht, the computa- 
tion of the likelihood function is mapped on to a similar 
problem for SV models, and hence existing techniques for 
further approximations are available. 

To our knowledge the present paper is the first to 
present results on ML estimation for multifractal mod- 
els with continuous state spaces for the volatility. Such 
estimates may be of great practical importance, since 
accurate parameter estimation is essential for volatility 
forecasts and risk estimates. In the MRW model this 
degree of accuracy is particularly important for the in- 
termittency parameter A which determines the peaked- 
ness of the return distributions on all time scales. In 
applications other than finance, accurate estimates of A 
can be used as supplements to the empirical scaling func- 
tions, and thereby the ML estimator can provide a tool 
for quantifying multifractality in data. 

The paper is organized as follows: In section [IT] we 
briefly explain the construction of the continuous-time 
process X{t) for which the model given by equations ^ 
and ([3]) is a discretization. There exists a large class of 
multifractal processes which are related to a construction 



known as infinitely divisible cascades (IDC). In general 
the random walk models associated with IDC processes 
have logartithmic volatility with infinitely divisible dis- 
tributions, and the MRW model considered in this paper 
is a discrete approximation to the random walk model ob- 
tained in the special case when the logarithmic volatility 
is Gaussian. 

Section |III| contains the procedure for approximated 
ML estimation in the MRW model. In section ITVl we test 
the estimator by first applying it to various stock market 
indices, and then by running a small Monte Carlo study. 
The results are compared with the GMM method used 
in p^. 

We finally remark that the methods presented in this 
paper have been implemented in a package for the R 
statistical software [20j. This package is available online 

m- 

II. MOTIVATION OF THE MODEL 

There exist several popular models for multifractal 
stochastic processes with uncorrelated increments. All 
of these models can be written either on the form 
X(t) ~ B(A{t)), where B{t) is a Brownian motion and 
A(t) = m{[0,t\) is the distribution function of a multi- 
fractal random measure m on the time axis, or as 

X{t) = lim / ^Ar(t')dB{t') , 
■•■^o Jo 

where Ar(t) A{t) as r — > 0. The meaning of Ar{t) 
is discussed below. These two types of constructions are 
equivalent as long as B{t) is a Brownian motion. (This is 
not the case for fractional Brownian motions with H 7^ 
1/2.) 

The differences between the various multifractal mod- 
els are then related to the construction of the random 
measure m. The log-normal MRW model is on one hand 
based on a particular construction of m known as multi- 
plicative chaos, and on the other hand it can be seen as 
a special case of the more general IDC constructions. 

In multiplicative chaos, which was first developed rig- 
orously in [3], one considers a sequence m„ of measures 
defined via random densities on the form 

dmn{t) = Cn e''"(*^dt, 

where l/c„ = E[e''"(*)], and /i„(i) are centered Gaus- 
sian processes with co- variance structures gn(t,s) = 
Cov{hn{t),hn{s)) that converge to some expression 
g{t, s) in the limit n — >■ 0. Kahane [3] showed that if g is 
CT-positive, meaning that there are positive and positive 
definite functions Km{t,s) such that 

n 

9n{t,s) = ^ K,n{t,s) , 
m— 1 

then the sequence to„ converges weakly to a Borel mea- 
sure m which depends only on the function g(t, s). One 
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can therefore informally think of m as being on the form 
dm{t) = ce'^^^Ut where 1/c = E[e''(*)] and h{t) is a 
"Gaussian" process with co- variance g{t, s). Then, if one 
makes the choice 

7(t,s) = A2log+-^, (4) 
\t-s\ 

one easily obtains the relation h{at) = h{t) + Cl{c), where 
n{a) are independent of h{t) and distributed according 

to ri(a) ^ A/'(0, -A^ logo). It follows that we iov t < R 
and < a < 1 have the scaling relation 

m{[0,at]) ^ M{a)m{[0,t]), (5) 

with log M(a) - A/'((l + AV2)loga,-A2loga). See 
proposition 3.3 in for a rigorous proof of and see 
example 2.3 of the same paper for a verification that the 
function g{t, s) in equation Q is cr-positive. By using the 
well-known formula for the q-th moments of log-normal 
variables together with equation ^ , we easily verify the 
multifractality of the process A{t) — m{[0,t]): Denote 
Cq = E|m([0, and observe that 

E|m([0,i])|« = CgEM{ty = Cgt^^'-'^^ , 

where (Aiq) =</(! + A^/2) — A^g^/2. Since a Brownian 
motion is self-similar with H = 1/2 the scaling function 
of X{t) is given by C(g) - a(g/2). 

Alternatively the model defined by equations ^ and 
([3]) can be motivated by considering the more general 
class of IDC models. Here we briefly mention the main 
ideas and results in this theory, and we refer to [23J HI] 
for details. At the base of this construction is an ob- 
ject called an independently scattered infinitely divisi- 
ble random measure P(dt, dr) defined on the halfplane 
S+ = {(i,r) e 1^ > 0}. The defining properties of 
the random measure P are: (1) for any measurable set 
A C 5*"*", the random variable P{A) is infinitely divisible 
with characteristic function 

where fj.{dt,dr) = r~'^dtdr. (2) for any finite sequence 
Ak C of disjoint and measurable sets, the correspond- 
ing random variables P{Ak) are independent. If we as- 
sume that ip'{0) — 0, the random measure P induces 
a family of centered and stationary stochastic processes 
through the equation 

Kit) = P{A{r,t)) , 

where A{r, t) are cone-likc domains defined by 

A{r,t) = {{t'y) e 5+ |r' > r, \t' - t\ < /(r')/2} , 

with f{r) = r ioi r < R and /(r) = R ior r > R. The 
time correlations in the processes hr(t) are characterized 



by the functions 

Prit) = M(^(0,r)n^(t,r)) 

^ flogf + , t<r 

nog+ f , t>r ' 

In fact, the co- variance of hr{t) is given by 

C0Y{hr{t),hr{s)) = X'prilt - S\) , 

where A^ = -^"{0). 

Random measures are defined by dmr{t) = c^e''''^*-' dt, 
where l/c^ — E[e'''^'-*^]. The corresponding distribution 
functions are Ar{t) ~ TOr([0,i]) and corresponding ran- 
dom walks are 

Xr{t) = [ ^/Ar{t')dB{t') . 

Jq 

By using the relation par (at) = — log a + pr (t) one can 
show that 

har{a-t) = hr{t) + Q{a) , 

for a e (0, 1) and t < R, where n{a) are independent 
of hr(t) and have characteristic functions ipn{a){Q) = 
g-V(9)ioga_ Consequently the limit process X{t) = 
Xr(t) has scaling function 

C(g)- (l + VH))'z/2-^(-*g/2). 

In the case that hr{t) are Gaussian, i.e. ip{q) = —X^q^/2, 
the co-variance is on the form 

for |i — s| > r, and hence it can be approximated by the 
process defined by equation In this case the scaling 
function is 

C(g) = (l + AV2) g/2-AV/8. 

We note that for A = the process X{t) is reduced to a 
Brownian motion and C.{q) = q/2. 

We point out that this paper presents a ML estima- 
tor for the discrete-time process Xt defined by equations 
([2| and (|3|. This is sufficient for the purpose of mod- 
eling and forecasting volatility in financial time series, 
since the discrete-time MRW model is directly compara- 
ble to GARCH-type models. In other applications, such 
as modeling the velocity field in turbulence, one is inter- 
ested in the continuous-time process X{t). Since xt is an 
approximation to the continuous-time process X(t), our 
method can also be interpreted as an estimator for this 
process. In this case one must be aware that the incre- 
ment process X{t+At)—X{t) is not proportional (in law) 
to e^^'*^*^et, and that this is only an approximation in the 
limit A^ <C 1. See appendix A.l. in [TOl. In the case of 
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strong intermittency, the estimator for the continuous- 
time process X{t) may therefore depend significantly on 
the time-scale At for which the data is sampled. An 
analysis of how our method performs as an estimator 
for X{t) will require extensive Monte Carlo simulations 
(with varying A and At) , and this is beyond the scope of 
this paper. 

We also remark that it in some applications is rele- 
vant to estimate the parameters of the measure dm{t), 
for instance when modeling the energy dissipation fields 
in turbulence. In the discrete-time approximation this 
corresponds to the process e''* , where ht is described by 
equation ([3|. Since ht is Gaussian, this problem is much 
easier than the one considered in this paper. The ML es- 
timator for e''* can be constructed using standard meth- 
ods [25j and no approximations are required. 



III. APPROXIMATED MAXIMUM 
LIKELIHOOD 

In this section we explain our method of approximated 
maximum likelihood estimation. Let xt and ht be the 
processes defined by ^ and Denote x = (xi,. . . , Xn) 
and h = (hi, . . . , /i„). The first step is to write 



Pxix) 



Px,h{x., h) dh 



Px\h{x\h)ph{h)dh. (6) 



The first factor Px\h{x\h) in the integrand is computed 
by noting that, when conditioned on h, the variables 
xi, . . . ,x„ are independent and Gaussian. In fact, 



t=i 

= — nlog\/27rccr 



(7) 



xt 



2 2cr2ce''t 



For the second factor Ph{h) we use that ht is a centered 
Gaussian process with a specified co-variance structure 
Cov{ht,hs) = j{\t — s\). First we decompose the density 
into one-dimensional marginals; 



Ph 



{h) ^ PhAhl)Wph,\h^,^^^{ht\hl■.t-l) , 



(8) 



t=2 



where we have used the notation 



(hn, hn+l, ■ ■ 
[hn 1 hn—1 : ■ ■ 



hm) for m > n 
for n > m 



■ ; hui 



Denote by Ft the co-variance matrix of the vector hi-^, 
and let 7i:t — (7(1), . . . , 7(<)). The co-variance matrix 
can be written on the block form: 



Ft = 



7(0) 7l:t-] 
lT:t-l Tt-l 



By performing standard computations of conditional 
marginals in multivariable Gaussian distributions we de- 
duce that htlhi-t-i is a Gaussian with mean 



rrit = 7l:t-irf_\/l(J_i):i 



and variance 



= 7(0) - Iv.t-iTi-ill.t-i ■ 

As usual it is convenient to introduce vectors 0^*^ defined 
by (/)'-*-'Ft = 7i:t. This allows us to write the mean as 
nit — 0*'*^^''^(t_i)-i ^-iid the variance as St — 7(0) — 
0(t-i)-y^^_^, Then from equation i8h we have 



(9) 



log Ph (h) = -n log V2tt - ^ log St 



2Sf 



Combining equation ^ with equation ([t]) we get an ex- 
pression for the full likelihood: 



logPx,h{x,h) = -nlog(27r Veer) 



E 



ht 



2cT2ce''* 



J2^og St 

i = l 



(10) 



252 



We keep in mind that c depends on R and A through 
the relation 1/c = E[e''*] = i?^'/^^ 

Approximation 1: By comparing co- variances the pro- 
cess ht can be written as 



<f>t'^ht^,+ 



^ti^h,+wt, (11) 



where Wt are independent Gaussian variables with zero 
mean and variances equal to St- As approximations to 
ht we can consider processes obtained by truncating the 
sum in equation (111. We fix a parameter r G N , and 



for t > T we replace equation (11 1 with 

ht = Cl)[^ht-l +■■■+ 4>i"^ht-r + w[^^ , (12) 
(t) 

where wj. are independent Gaussian variables with 
zero mean and variances equal to S^^i- Note that 



ht\hj_i.t_^ in equation (11) has the same distribution 



as obtained from (12), namely a Gaussian with mean 
nit — '^'■^''^(t-i) t-r ^^"^ variance S'^^i- In effect we have 
approximated the distribution of ht\ht-i:i, by truncat- 
ing the dependency after a lag r. As a result of this 
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approximation equation ^ becomes 
\ogph{h) = -nlogV2TT 

r 

- ^ log St - {n - t) log Sr+i 
"-"I'f.-iyif 



(13) 



2Sf 



E 

t=T + l 



In order to compute the expression in equation ( 13 ) we 
need to solve the equations 



(*)r, = 



7l:t , i = 1, 



This is done efficiently using the Durbin-Levinson algo- 
rithm [26l [27], W e remark that for r = n the expression 
in equation (13) is exact. 



Approximation 2: The second approximation is the 
so-called Laplace's method, which is frequently used for 
approximation of likelihoods in SV models, see e.g. [551 
[55]. We write equation ^ on the form 



(14) 



where 



n 



1 

^^logp^^l^^(xt) 

t=i 
1 " 



(15) 



Laplace's method is to assume that fx{h) has a global 
maximum in M", which we denote by h* . When n is 



large the contribution to the integral in equation ( 14 ) is 
concentrated around h* , and therefore we make a second 
order Taylor approximation to fxih) around this point. 
Since h* is also a local maximum we have 

Uh) « - logp,.^(a;, h*) + ^{h- h*) (h - h* f , 
n 2n 



where 



^x — 



d'^ ^ogPx,h{x,h*) 
dhdh^ 



is the Hessian matrix of ,fx(h) at the point h* . The ap- 
proximation now reads 

= (2^)"/2|detf],|-i/2p,,;,(x,/i*). 

The maximum h* is found by computing the partial 
derivatives of fxih) with respect to h, setting them equal 
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FIG. 1. The top figure shows the daily log-returns of the Ger- 
man DAX index for the time period 1990/11/26-2011/11/25. 
The standard deviation of the data is normalized to unity. For 
T — 500 the ML estimates are A = 0.32 and T = 2.2 years. 
The lower figure shows a simulation of the MRW model Xt 
with the estimated parameters. 



to zero and solving the corresponding system of equations 
numerically using the algorithm DF-SANE [3D], which is 
implemented in R package "BB" [31]. The matrix 
is obtained using analytical expressions for the second 
derivatives. This matrix is band-diagonal with band- 
width equal to the truncation parameter r, and in the 
R software such matrices are efficiently stored and ma- 
nipulated using the package "Matrix" [32]. 



TABLE I. Estimated parameters for the log-returns of various 
stock market indices. Prior to the analysis the sample stan- 
dard deviation of each data set is normalized to unity. All 
ML estimates are run with r = 500 and the GMM estimates 
are performed with a maximum time lag tmax = 500 days in 
the auto-correlation function of nit = loga;^. The analyzed 
data is retrieved from 1331. 





ML 


GMM 


Index (time period) 


A 


T (years) 


A 


T (years) 


GAG 40 (1990-2011) 


0.29 


2.8 


0.36 


2.5 


S&P 500 (1950-2011) 


0.32 


12.2 


0.36 


10.2 


DAX (1990-2011) 


0.32 


3.3 


0.44 


3.1 


Nikkei 225 (1984-2011) 


0.36 


1.4 


0.40 


3.0 


Hang Seng (1986-2011) 


0.37 


2.8 


0.44 


2.5 


FTSE 100 (1984-2011) 


0.28 


4.2 


0.36 


2.9 



IV. ESTIMATOR COMPARISONS 

In this section the ML estimator is compared with an 
GMM approach which is similar to the one used in [19j . 
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PDF of GMM estimate for A. 
Sample length n=5000. 



Li. 



0.2 0.4 0.6 

PDF of GMM estimate for A. 
Sample length «= 10000. 
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PDF of GMM estimate for tr. 
Sample length n=2500. 
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PDF of GMM estimate for cr. 
Sample length n=5000. 
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PDF of GMM estimate for o". 
Sample length «=10000. 
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(c) PDF of GMM estimate for log R. 
Sample length n=2500. 
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(f) PDF of GMM estimate for log R. 
Sample length «=5000. 
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(i) PDF of GMM estimate for log R. 
Sample length n= 10000. 
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FIG. 2. The results of the Monte Carlo study for the GMM estimator explained in section [TVl The figures show the estimated 
probability density functions for the estimators based on 500 realizations of the process. The parameters are A — 0.35, cr = 1 
and R = 2000 (i.e. logR — 7.6). In figures (a-c) the sample lengths are n = 2500, in figures (d-f) the sample lengths are 
n = 5000 and in figures (g-i) the sample lengths are n = 10000. The means and standard deviations of the estimators are 
reported in table [III 



This GMM version is essentially a least-square fitting of 
the auto-correlation function for the logarithmic volatil- 
ity, and we briefly explain this method in the following: 
Denote rrit — log Xf and observe that 

■mt = ht+ yt 

where yt — log c + log £\ are independent and identically 
distributed. We can use the sample standard deviation 
to normalize rrit so that it has unit variance. Then, if 
we let /i„j = E[mt] = E[yt] denote the mean of to^, the 
auto-correlation function of rrit has the form 

ACF„j(t) = E[(toi - ^„)(mt+i - ^™)] 
= E[/ii/it+i] =A2log+-^. 

For t < R we have 

ACF™(i) = A^logi?- A2log(t + l), 

and log R and A can be found by linear regression of the 
auto-correlation function versus log{t -1-1). 

We begin testing the approximated ML estimator by 
applying it to various stock market indices. We use daily 



log-returns and in all of the estimates the truncation pa- 
rameter is set to T = 500 days. The results are presented 
in table |I] We observe that the intermittency parame- 
ter A varies from 0.29 to 0.37 for the different indices 
and time periods. We also observe that the correlation 
range parameter T varies by roughly one order of mag- 
nitude, in the range 1.4-12.2 years. If we compare with 
the GMM we see that, for all the indices, the estimates 
of A are lower using the ML method. For the parameter 
T the estimates using ML and GMM are more or less 
consistent, but with quite large variations between the 
two estimators. 

To further test the performance of the proposed ML 
estimator we run a small-sample Monte Carlo study. 
We have used three different sample lengths n S 
{2500,5000,10000}, and for each sample length n we 
simulated 500 sample realizations. The parameter vec- 
tor considered is A = 0.35, a — 1 and R = 2000. For 
the truncation parameter we have considered the cases 
T G {10, 50, 100}, and in the GMM method we use a max- 
imum time lag imax = 500 days in the auto-correlation 
function of rrit = log . 

The results are presented in table |llj For both the 
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(a) PDF of ML estimate for A with T= 100. 
Sample length «=2500. 
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(d) PDF of ML estimate for A with t=100. 
Sample length j!=5000. 
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FIG. 3. The results of the Monte Carlo study for the ML estimator with r = 100. The figures show the estimateii probability 
density functions for the estimators based on 500 realizations of the process. The parameters are A = 0.35, a = 1 and R — 2000 
(i.e. logi? — 7.6). In figures (a-c) the sample lengths are n — 2500, in figures (d-f) the sample lengths are n = 5000 and in 
figures (g-i) the sample lengths are n = 10000. The means and standard deviations of the estimators are reported in table [ill 
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FIG. 4. Double-logarithmic plot of the mean square errors 
as functions of sample length n for the ML estimator with 
r = 100 (squares) and the GMM estimator (crosses). The 
dotted lines have slopes equal to —1, i.e. the mean square 
errors decay roughly as 1/n for both estimators. 



GMM and the ML methods the estimates of R are highly 
unstable. This is also pointed out in [TH]. However, the 
processes Xt only depend on the R through expressions 
on the form log R. Therefore, in order to have an esti- 
mator which is comparable to A, we should consider the 
variable logi?. The estimators of logi? behave reason- 
ably well, even though there are significant mean square 
errors and some bias. We see that both the ML and 
GMM method underestimate logi? and that the errors 
are roughly the same for the two estimators. 

On the other hand we observe that the ML estimates 
of A have a standard deviations which are much smaller 
than the corresponding standard deviation for the GMM 
estimate, especially for r = 100. This can also be seen 
from figures [2] and [3] where the probability density func- 
tions for the GMM estimates and the ML estimates are 
presented. Based on this we conclude that the ML esti- 
mator performs better than the GMM. Moreover, if one 
allows longer computing times, the truncation parame- 
ter T can be increased to obtain even more accurate es- 
timates. For a time series oi n — data points, an ML 
estimate with r = 500 takes a few minutes on a personal 
computer. 

In figure [4] we have plotted the mean square errors 
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TABLE IL The results of a Monte Carlo study of the ML and 
GMM estimators. The parameters in the simulations are A = 
0.35 and R = 2000 (i.e. logT? = 7.6). In the GMM estimator 
we have used a maximum time lag tmax = 500 days in the 
auto-correlation function of mt = loga:^ . The reported values 
are the mean estimates together with the standard deviations 
(in brackets). 





ML 


GMM 


n 


T 


A 


logi? 


a 


A 


logi? 






10 


0.31 
(0.03) 


6.87 

(3 4Ll 


0.97 

("0 19") 








2500 


50 


0.34 
(0.03) 


6.47 
(1.73) 


0.97 
(0.19) 


0.34 
(0.08) 


6.11 
(0.76) 


0.97 
(0.19) 




100 


0.34 
(0.03) 


6.35 
(1.67) 


0.97 
(0.19) 










10 


0.30 
(0.03) 


5.58 
(2.18) 


0.98 
(0.14) 








5000 


50 


0.34 
(0.02) 


7.02 
(1.44) 


0.98 
(0.14) 


0.35 
(0.05) 


6.69 
(0.96) 


0.981 
(0.15) 




100 


0.34 
(0.02) 


6.87 
(1.31) 


0.97 
(0.14) 










10 


0.30 
(0.02) 


9.10 
(1.80) 


0.98 
(0.10) 








10000 


50 


0.34 
(0.01) 


7.37 
(1.24) 


0.98 
(0.10) 


0.35 
(0.04) 


7.11 
(0.92) 


0.98 
(0.10) 




100 


0.34 
(0.01) 


7.21 
(1.16) 


0.98 
(0.10) 









(MSE) E[(A - \f] for the ML estimator with t = 100 
and the GMM estimator. We see that for both the esti- 
mators the MSE is roughly invserly proportional to the 
sample length. However, from table |lT] we see that is a 
slight negative bias in A for the ML estimator. This bias 
decreases with increasing r, and we suspect the estimator 
to be asymptotically unbiased in the limit t ~ n ^ oo. 



V. CONCLUDING REMARKS 

In this paper we have presented an approximate ML es- 
timator for MRW processes. The method is implemented 
and tested in a Monte Carlo study, and the results show 
significant improvements over existing methods for the 
intermittency parameter A. 

The methods of this paper represent a suitable start- 
ing point for two important generalizations. The first 
generalization is to allow for correlated innovations, for 
instance by letting et be a fractional Gaussian noise. This 
has several important applications, for instance in mod- 
eling of geomagnetic activity [S^ rf and electricity spot 
prices [331 . Another interesting generalization is to con- 
sider the non-Gaussian IDC models referred to in section 

im 

We also point out that the techniques presented in 
section Hill can be used to calculate conditional densities 
on the form p{xt+i, ■ ■ ■ ,xt+s\xi, . . . ,Xt)- At time t, 



such an expression provides a complete forecast over the 
next s time steps. Forecasting and risk analysis based 
on the MRW model and the methods in this paper is a 
promising topic that will be pursued in future work. 
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